// Special version of run_QN to test pair_potential_function and
// allow a deformation process of this function

#include <iostream>
#include <fstream>
#include <sstream>
#include <cmath>
#include <cstdio>
#include <string>
#include <iomanip>

#include "pair_potential_function.h"

#include "continuous_tabusearch.h"
#include "Time.h"
#include "wrapperFunctions.h"
#include "wrapper_quasi.h"

#define PI   3.1415926535897932385

using namespace std;

struct StartEntry
{
  double value;        // fom value 
  vector<double> coor; // coordinates
};

// to read in starting solutions 
void readStartingSolutions(vector<StartEntry>& data)
{
  ifstream in("./input/startingConfigurations5.txt");

  if ( !in.is_open() )
  {
    cout << "can't open file ./input/startingConfigurations.txt \n";
    exit(1);
  }

  if ( data.size() != 0)
  {
    cout << "\nIn run_QN_time.cpp assume data is empty\n";
    exit(1);
  }

  string aLine;
  int line = 0;
  while ( getline(in, aLine, '\n') ) 
  {
    line++;

    stringstream str(aLine);
    double value;

    StartEntry entry;
    int index = 0;
    while ( str >> value )
    {
      if (index == 0)
      {
        entry.value = value;
      }
      else
      {
        entry.coor.push_back(value);
      }
      index++;
    }
    data.push_back(entry);   
  }
}

void run_QN_time(int argc, char *argv[ ], int seed, double gtol, Function *func, int nruns, int numdim, string& directory)
{
  ofstream fpParams("./output/QN_params.used");
  fpParams << "Seed               " << seed << endl;
  fpParams << "Tolerance          " << gtol << endl;
  fpParams << "N_runs             " << nruns << endl;
  fpParams << "N                  " << numdim << endl;

  fpParams.close();

  initiate_ran1(seed);      // initiates the random number

  Time timer;               // to measure calculation times
  cout << "prepare timer " << timer;

  double *sol;              // the solutions we will follow
  //double *deriv;
  double gTol = gtol;
  int nRuns = nruns;
  int numDim = numdim;


  ////////////////////////////////////////////////////////////////
  // Open files QN_Solution.out and QN_initial_values.out,
  // and prepare them for output. QN_Solution.out holds the
  // best function value after each run both without local
  // minimisation and with local minimisation. QN_initial_values
  // holds the initial coordinates (debugging).
  ////////////////////////////////////////////////////////////////

  ofstream outEndQN_Solution( string(directory + "QN_Solution.out").c_str());

  if ( !outEndQN_Solution.is_open() )
  {
    cout << "Cannot open file in run_QN.cpp\n";
    exit(1);
  }

  outEndQN_Solution << "Dimension: " << numDim << endl
    << "Smallest function value after each run\n\n"
    << "      Func Val Bef     Func Val After   Function Evaluations      Distance       Run Number" << endl
    << endl;

  // output final values in one file which can easily be read in 
  // to form the starting values from which to start deformation 
  // approach

  ofstream outFinalVal("output/QN_final_values.out");

  // other files to output

  ofstream outIniVal("output/QN_initial_values.out");
  ofstream outIniValCif("output/QN_initial_CIF.cif");
  ofstream outFinalValCif("output/QN_final_CIF.cif");


  Function *function;       // Creates an instance of the function
  function = func;          // Sets the function = to the function passed as an argument

  Funcdd funcd(numDim, function);

  newA(sol, numDim);

  vector<StartEntry> data;
  readStartingSolutions(data);
  if (numDim != data[0].coor.size())
  {
    cout << "\nDimension mismatch. numDim = " << numDim << " data coor vector size = " 
      << data[0].coor.size() << endl;
    exit(1);
  }

  cout << "Time: " << timer;

  PairPotentialFunction* ppf = dynamic_cast<PairPotentialFunction*>(function);
  ppf->setDeformation(5);

  for(int i=1; i<=nRuns; i++)  
  //for (unsigned int i=0; i<data.size(); i++) 
  //for (unsigned int i=0; i<3; i++)
  {
    cout << "Run No. " << i << endl;

    function->nCall = 0;
    // Populate the Input array with a random start position/configuration
    function->createRandomSol(sol);
    //for (int j = 0; j < numDim; j++)
    //  sol[j+1] = data[i].coor[j];


    // Print out the run number and initial
    // starting configuration (coordinates)
    outIniVal << "\nRun No. " << i << endl << endl;
    outIniVal << "Variable | Starting values" << endl;
    for (int k=1; k<=numDim; k++)
    {
      outIniVal << setw(5) << k << setw(18) << sol[k] << endl;
    }           

    double origValue = function->Value(sol);

    // Also write out initial starting configs as cifs
    stringstream label_ini;
    label_ini << i << "_" << origValue;
    function->print_sol(sol, origValue, outIniValCif, string("Run_") + label_ini.str());
    outIniValCif << endl;

    // print out best solution so far (without local minimisation)
    outEndQN_Solution << setw(16) << origValue;

    

    // output follow path
   /* stringstream str;
    str << directory << "followPath_" <<i << ".cif";
    ofstream followPath( str.str().c_str() );

    // do deformation
    for (double t = 0; t <= 400; t+=0.1)
    {
      cout << "t = " << t << endl;
      ppf->setDeformation(t);
      origValue = function->Value(sol); */
      WrapperQuasi(sol, gTol, origValue, funcd, directory); 
/*
      // output follow-path
      stringstream label_path;
      label_path << t << "_" << origValue;
      function->print_sol(sol, origValue, followPath, string("Entry_") + label_path.str());
      followPath << endl;
    }
*/
    // note after local minimisation origValue has been updated
    outEndQN_Solution << setw(18) << origValue << setw(23) << function->nCall << setw(19)
      << function->distanceFromGlobalMinimum(sol) << setw(14) << i << setw(14) << gtol << endl;

    // Also write out final configs as cifs
    stringstream label_final;
    label_final << i << "_" << origValue;
    function->print_sol(sol, origValue, outFinalValCif, string("Run_") + label_final.str());
    outFinalValCif << endl;

    // write out file for created starting values from which to do deformation from
    outFinalVal.precision(5);
    outFinalVal << setw(17) << origValue;
    outFinalVal.precision(9);
    for (int i = 1; i <= numDim; i++)
      outFinalVal << setw(17) << sol[i];
    outFinalVal << endl;
  }

  double *bestSol;
  newA(bestSol, numDim);
  char resultBest[80];

  sprintf_s(resultBest, 80, string(directory + "QN_best_struc").c_str());
  double bestVal = function->BestSoFar(bestSol);
  //function->print_sol(bestSol, bestVal, resultBest, 0, 0);

  outEndQN_Solution << "\nGlobal best value over all runs:\t" << setw(9) << function->BestSoFar() << endl;

  outEndQN_Solution.close();
  //outIniVal.close();

  // print final output to screen
  cout << endl << endl << "Output files have been written to " << directory << endl;
  cout << "The total time for the simulation was " << timer << endl << endl;
  cout << "---------------------------" << endl << endl;
}

